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Abstract 

We present upper bounds for the numerical errors introduced when using operator splitting methods to integrate 
transport and non-linear chemistry processes in global chemical transport models (CTM). We show that (a) operator 
splitting strategies that evaluate the stiff non-linear chemistry operator at the end of the time step are more accurate, 
and (b) the results of numerical simulations that use different operator splitting strategies differ by at most 10%, in a 
prototype one-dimensional non-linear chemistry-transport model. We find similar upper bounds in operator splitting 
numerical errors in global CTM simulations. 
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1. Introduction 


Global tropospheric chemistry transport models 
(CTM) are used to address important issues ranging 
from air quality to climate change. In order to con¬ 
tinuously improve their performance, it is of crucial 
importance to understand and quantify the diverse 
sources of uncertainties and errors present in them. 
We group these in three different categories, (i) errors 
and uncertainties coming from observations and data 
used in our models (such as emission inventories, 
wind fields, reaction rates); (ii) errors coming from 
our choice of governing equations (or mathematical 
model), parametrizations, and the level of complexity 
of the physical modules included in our formulation; 
and (Hi) numerical errors coming from the choice of 
algorithms we use to solve the governing equations 
using computers (Enting, 2002; Zhang et al, 2011). 


In this study, we focus our attention on estimating the 
magnitude of numerical errors (Hi), in particular, those 
arising from the choice of operator splitting technique 
utilized to integrate in time the transport and chemistry 
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operators in real-life global CTMs. In order to achieve 
this, we numerically extend the results introduced for 
the linear diffusion-reaction case in (Sportisse 2000| ), 
to a non-linear 1-D chemistry-transport numerical 
model. The latter numerical results provide us with a 
framework to estimate upper bounds for operator split¬ 
ting errors in the fully non-linear 3-D state-of-the-art 
global CTM: GEOS-Chem ( ]Bey et al.[ (200T ). To the 
best of our knowledge, our contribution is the first in 
estimating operator splitting errors in the context of 
real-life global atmospheric chemistry simulations. 


CTMs simulate the dynamics of chemical species in 
the atmosphere by numerically integrating a set of cou¬ 
pled nonlinear partial differential equations of the type: 


^+V-(h Cd = V-(pKV < nyP i (C j )-C i L i (Cj)+Qi-S i 

( 1 ) 

for i = 1,...,A; where Ci(x,t) represents the spatio- 
temporal evolution of the concentration of species 
i (typically over a hundred species are considered), 
u(x , t) is the wind velocity, p is the air density, K the 
eddy diffusivity matrix, P t are the nonlinear production 
terms, L, are the destruction terms, Q t are the volume 
emission sources, and Si are the sinks (ex. precipitation 
or in-cloud removal). See |Sportisse| ( J20Q7| ) for a 
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detailed description of these equations. 


Due to the dimensions of grid boxes in global CTMs, 
like GEOS-Chem (with hundreds of kilometers in the 
horizontal versus tens to hundreds of meters in the 
vertical), intertial vertical transport processes in this 
global models are simulated (a) using vertical mass 
fluxes schemes that ensure that the horizontal air flow 
is divergent-free (V/* or u = 0), (b) using convection 
parametrizations, and (c) using a boundary layer mix- 
ing algorithm (|Lin and Rood 1996 Allen et ak||1996 


|Wild and Prather[ 2006HPrather et al.' |2008 ). In ad 

dition, horizontal diffusion due to numerical errors in 
transport schemes are typically higher than their Eddy 
difusivity counterpart, as measured by aircraft missions 


([Pisso et al. | 

2009 [ 

Wild and Prather, 2006; Rastigeyev 

et al. 2007; 

Santill 

ana 2013). As a consequence, the 


models the dynamics of intertial vertical transport as 
an eddy diffusion process, is not explicitly integrated 
in global CTMs; and the governing equati ons ([T]) are 


sometimes written ( Rastigeyev et al.| [2007 1 |S antiliana 
et al.| |2010[ |Santillana 2013 [ ) in a simplified way as 


+ m • VC* = Pi(Cj) - CiU(Cj ) + Qi - Si. (2) 


The chemistry operator on the right-hand-side of 
equations 0 models the chemical interaction of 
atmospheric species whose lifetimes range from mil¬ 
liseconds to many years. The chemistry operator is very 
stiff as a consequence of this large range of time-scales 
and thus, implicit-in-time methods are an appropriate 
choice to integrate equations 0 - Traditional methods, 
such as the method of lines, aimed at achieving this 
task in realistic 3D simulations, involve solving for an 
enormous number of degrees of freedom at each time 
step in a coupled fashion ( 10 8 « 100 chemical species 
in ~ 10 6 grid cells, for a 1° x 1° spatial resolution). This 
is due to the inter-species coupling in the chemistry op¬ 
erator and the spatial coupling in the transport operator. 
In practical situations, however, efficient computational 
algorithms to integrate equations 0 use operator split¬ 
ting strategies that allow the explicit time-integration 
of the transport and implicit time-integration of the 
chemistry operators separately and sequentially, thus, 
reducing significantly the degrees of freedom solved in 
a coupled fashion at a given time step. This is done at 
the expense of a loss of accuracy in the approximate 
solution ( [Hundsdorfer and Verwer||2003| . 


in realistic 3-D computer simulations is a hard task 
since no relevant analytic solution can be used as a 
reference to estimate them. In theory, estimates of these 
errors depend directly on the regularity properties of 
the analytic solution of equations ([T}, the set of initial 
and boundary conditions, and the chosen numerical 
scheme ( |Guo and Babuska"| |1986[ |Iserles[ |2QQ9| |Ern| 
land Guermond| [20041 |Brenner and Scott||2Q08| ). In this 
study, we assume that the analytic solution of equations 
0 is unique and regular enough so that numerical 
error estimates can be expressed as inequalities of 
the form 0 - Operator splitting errors, as well as 
numerical errors arising from the time-integration 
of the chemistry operator depend explicitly on the 
magnitude of the chosen time steps, while numerical 
errors coming from the time-integration of the transport 
operator depend both on the time step and on the grid 
size. This fact, in combination with an expression of 
the analytic solution of equations 0 . is exploited to 
obtain the exact magnitude of operator splitting errors 
in our one-dimensional proto-type transport-chemistry 
numerical model. 


Our one-dimensional numerical experiments show 
three main results: (a) operator splitting sequences 
where the stiff non-linear chemistry operator is evalu¬ 
ated at the end of the time step are more accurate than 
those where the transport is evaluated lastly, indepen¬ 
dently of the operator splitting time-step, as in the linear 
case introduced in ( |Sportisse[ |2000| ); (b) the results of 
numerical simulations that use different operator split¬ 
ting strategies differ by at most 10%; and (c) numer¬ 
ical errors coming from the integration of the trans¬ 
port operator are much bigger than those coming form 
the operator splitting technique for spatial and temporal 
scales comparable to those used in global CTM. We use 
this fact, and evidence from papers suc h as ([Wild and 


2008 


Santillana 


Prather 2006[ Rastigey ev et al.| |2007| |Prather et al 


2013| , to suggest that in realistic 3D 


simulations, errors due to operator splitting are much 
smaller than those introduced by transport schemes. 


2. Numerical error estimation 

Upper bounds of the numerical errors introduced 
by solving partial differential equations with regular 
boundary and initial conditions, using a given numerical 
scheme, can be expressed by inequalities represented as 

|| C(x, t) - C h (x, 0ll Vl < Mi At a + M 2 Ax p (3) 


Estimating the magnitude of the numerical errors 
introduced by the time-integration of equations 0 


where C(x, t) is the true solution of the partial differ¬ 
ential equation, Ch(x,t) the numerical approximation, 
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At and Ax are the time step and grid size respectively, 
a and y 3 are exponents (typically larger than one) that 
determine the order of convergence of the method in 
time and space respectively, M\ and M 2 are constants 
that depend on the regularity of the true solution C(x, t) 
and parameters in the equation, and || • \\ V] is the norm 
in the appropriate Banach space V\. For a convergent 
method, as At —> 0 and Ax —> 0, the numerical error 
vanishes, (i.e. \\C - Ch\\ v —> 0) and the numerical 

approximation Ch converges to the true solution C, in 
the normed space V\. More details about the integral 
representation (equation [3} of numerical errors due to 
discretization of partial differential equations can be 
found in: Guo and Bab uska ]( [1986T ); Iserles|( [2009] ); |Ern| 
|and Guermond| ( [20Q4| ); |Brenner and Scott] ( |2QQ8| ) 

For the specific set of partial differential equations 
Q> operator splitting errors and errors coming from the 
numerical integration of the chemistry operator (where 
no coupling in space exists) contribute to the first 
term on the right-hand-side of inequality ([3]), whereas, 
numerical errors from the integration of the transport 
operator contribute to the first and second terms of 
the right-hand-side of inequality 0 - Quantifying the 
independent contribution of each processes to each term 
of inequality ([3} is not simple in practical applications. 
In the following section, we show how to estimate the 
magnitude of operator splitting errors in the absence of 
other numerical error coming from the time-integration 
of the transport and chemistry operators. 


2.7. Operator splitting techniques and error estimation 

Classical approaches to estimate the numerical errors 
introduced by operator splitting approaches are based 
on asymptotic expansions of exponential operators (lin¬ 
ear case) and Lie operator formalism (nonlinear case). 
For completeness, we briefly describe important results 
of the linear analysis of operator splitting methods in 
this section. We refer the reader to ILansen and Verwerl 
1999); |Sportisse| ( |2000| ); |Hundsdorfer and Verwer 


(2003) and the references therein for more details. In 


this section, it is assumed that the time-integration of 
each operator separately can be found exactly giving 
rise to no numerical error, i.e. the numerical errors 
discussed below come only from the choice of the 
operator splitting technique. 


We use as an example the linear evolution equation, 


where A and B are linear operators. One of these op¬ 
erators could represent the linear spatial differential op¬ 
erator d/dx (transport) in equations 0. The analytic 
solution for this problem is given by: 

v = vo exp((A + B)t) (5) 

The simplest operator splitting method, called Godunov 
and denoted by (A - B), can be obtained for t e [0, At] 
by solving the two evolution equations in sequence as: 

= Av*, v*(0) =vo in [0, At] 

dt 

dv** 

—— = £v**, v**(0) =v*(A t) in [0, At], 

dt 

( 6 ) 

The value for v at t = At is given by VAB(At) = 
v**(A t). The solution obtained with this operator split¬ 
ting method at t = At is given by 

VAB(At) = vo exp(BAt) exp(AAf) (7) 

The exact solution 0 and the solution vab in the previ¬ 
ous equation will be the same if 

exp((A + B)At) = exp(7?At)exp(AA0. 


This will happen if the operators A and B commute 
(think of matrices), i.e. if AB = BA. When AB ± BA , 
then the (point-wise) local-in-time numerical error asso¬ 
ciated to solving problem 0 using Godunov’s operator 
splitting technique can be shown to be 


leAB - 


(AB - BA) 2 
- z -Ar 2 v 0 


( 8 ) 


which leads to a global error 0(At ), i.e. 
||v — vab\\ < Mab At (for a constant Mab that de¬ 
pends only on the regularity of the analytic solution 
v). Since the numerical error vanishes as At —> 0, 
Godunov’s method is a convergent first order method 
in time, in the linear case. Another simple Godunov 
operator splitting can be obtained by reversing the 
order of evaluation of the operators A and B to obtain 
the (B - A) method ( v B a )• A more accurate and 
symmetric operator splitting method, often referred to 
as Strang method (Strang, 1968), can be obtained by 
averaging the output of the two previous methods, i.e. 
vs (At) = \(vab + vba)- It can be shown that the Strang 
operator method is globally second order accurate, i.e. 

11v — vs || < Ms At 2 for a constant Ms (Sportisse 


[Hundsdorfer and Verwer[|2003| ). 


2000 


The linear analysis presented above may fail and lead 
to different convergence results if one of the operators 


— = Av + Bv , v(0) = vo, 

dt 


V e IT (4) 
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is stiff, i.e. if the dynamics of one operator take place in 
much faster time scales than the dynamics in the other 
operator ( [Sportisse 2000| ). This can be seen by intro¬ 
ducing a small parameter e (representing the ratio be¬ 
tween fast time scales in the stiff operator and the slow 
time scales of the other operator) and re-writing the 
linear evolution equation ([4]) as a singular perturbation 
equation by re-defining 


use of implicit schemes to integrate the chemistry 
operator in global chemistry models leads to the choice 
of large operator splitting time-steps compared to the 
intrinsic stiffness of the chemistry system (At » e). As 
a consequence, and according to expression ( fT0| ), we 
may expect to observe large operator splitting errors 
when solving equations 0 with stiff and potentially 
non-linear chemistry operators. 


, x (<0 
6 


and 


B = T. 


(9) 


For our purposes, one can identify the chemistry opera¬ 
tor with the stiff operator^/6, (the nonlinear chemistry 
can be, locally-in-time and space, approximated by a 
linear and stiff mechanism at least for some subset of 
fast species), and identify the transport operator with 
the slow operator T, for which the dynamics takes place 
in a more confined range of time scales (as represented 


by our global models). It is shown in (Sportisse 2000) 
that the local error for the - T^j Godunov method be¬ 
comes (compare to equation ([8|): 


le £ 


(XT-T X ) 


Arv 0 


( 10 ) 


leading to a global error O(^), implying that 

||v - v 6 || < M 6 (y)- Note that convergence of the 
operator splitting method, in this case, can only be 
guaranteed provided the operator splitting time step, 
At, is small enough to satisfy At <$c e so that higher 
order terms, 0(^) k , will indeed vanish as k —> oo in the 
Taylor expansion of the error. 

In atmospheric chemistry simulations, we use oper¬ 
ator splitting methods to integrate in-time two opera- 


It is argued in [Sportisse]|2000| , that operator splitting 
errors, even in the presence of large operator splitting 
time steps (such that At » e ), may not be as big as sug¬ 
gested by expression |Sportisse| ( |2000| ) argues that 
the stiffness of the system can be balanced by the exis¬ 
tence of an underlying reduced model (low-dimensional 
manifold) describing the dynamics of the system and 
thus, by choosing the appropriate order of operator 
evaluation in a time-step, the splitting error may be 
bounded even with the increase of stiffness. Moreover, 
he shows for the linear case that sequences where 
the stiff operator is evaluated at the end of the time 
step lead to convergent and accurate methods in a one 
dimensional diffusion-chemistry toy example, even for 
large operator splitting time steps. In solving equations 
0, examples of these sequences include: Transport- 
Chemistry and Chemistry-Transport-Chemistry. 

Intuitively speaking, evaluating the transport operator 
at the end of the time step sets the state of the system 
far from the underlying low dimensional manifold 
driving the chemical system and provides an initial 
condition vo for the next time evaluation that enhances 
error propagation. This is avoided by evaluating the 
stiff chemistry operator at the end of the time step. The 
existence of these reduced models driving the dynamics 


tors in equations ([!]): transport and chemistry. Trans¬ 

in regional and global atmospheric chemistry models 

port and chemistry are known to commute when the ve¬ 

has been found in Lowe and Tomlin ( 

2000); Santillana 

locity field is divergent free and chemistry is indepen¬ 

et al. (2010|; Rastigeyev et al. (2007 

), suggesting that 

dent of the spatial location. In real atmospheric situ- 

the operator splitting order should be selected carefully. 


ations, these conditions are typically not met. Indeed, 
the non-linear chemistry operator depends dynamically 
on the geographic location (due to photolysis), and at¬ 
mospheric wind fields are in general not divergent-free. 

The result of the linear analysis above suggests that 
operator splitting approaches will converge only if the 
operator splitting time step is much smaller than the 
lifetime of the fastest species in the chemistry mech¬ 
anism (At 6). This is also the criterion established 
to ensure stability and convergence of explicit-in-time 
chemistry solvers, and suggests the use of prohibitively 
small operator splitting time-steps in order to guarantee 
convergence of the method. In practice, however, the 


To the best of our knowledge, a careful investigation 
of these errors in the realistic non-linear case does not 
exist so far and thus we aim at achieving this here. 

Isolating operator splitting errors in practical global 
atmospheric chemistry models is not straightforward, 
first, because we lack expressions for the analytic solu¬ 
tion of the system in realistic circumstances, and sec¬ 
ond, since the solutions of the chemistry and trans¬ 
port operators, separately, are obtained using numeri¬ 
cal schemes and thus are not exact as it was assumed in 
the previous analysis. In order to estimate upper bound 
estimates of operator splitting errors we proceeded as 
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follows. We first found sharp estimates of numerical er¬ 
rors in a ID non-linear chemistry-transport prototype 
problem with a known analytic solution. We designed 
this ID problem to resemble the interaction of numer¬ 
ical errors in the time-integration of the transport and 
(stiff) non-linear chemistry, when using operator split¬ 
ting methods, at spatial and time scales used in 3D 
global simulations. Our ID findings guide our method¬ 
ology to understand the differences observed between 
the outputs of 3D global simulations using different op¬ 
erator splitting strategies. We performed multiple 3D 
global simulations in order to further understand addi¬ 
tional numerical errors, due to the time integration of 
relevant processes (emissions, convective transport, and 
deposition) inherently solved with operator splitting ap¬ 
proaches. 


3. One-dimensional advection-reaction system 

We considered a one-dimensional advection-reaction 
system that can be solved analytically and thus exact 
values of numerical errors can be obtained. The system 
is characterized by a constant wind field throughout the 
domain, and a three-species ( NO , NO 2 , O 3 ) stiff non¬ 
linear chemistry-mechanism modeling the NO x (NO + 
NO 2 ) cycle through oxidation by ozone (0 3 ). This cycle 
is key in determining the balance of Ozone (0 3 ) in the 
atmosphere. The chemical reactions are given by: 

NO + 0 3 N0 2 , N0 2 ^ NO + 0 3 (11) 


where the parameters k\ and k 3 represent the constant 
reaction rates throughout the domain. The resulting 
advection-reaction system of equations can be written 
as 

d NO d NO 

—«-1" M — - = 

dt dx 

d N0 2 d N0 2 

dt U dx 

d0 3 d0 3 

dt U dx 

where NO, NO 2 , and O 3 , represent the concentration 
of each chemical in space and time, and u the constant 
velocity of the flow (compare with equations ([!])). 


-ki(NO) 0 3 + k 2 N0 2 (12) 

= k x (N0)0 3 -k 2 N0 2 (13) 

-ki(NO) 0 3 + k 2 N0 2 (14) 


The advection and reaction operators commute in 
this problem (since the advection operator is divergent- 
free, du/dx = 0, and the chemistry is independent of 
the location in space), thus, the use of operator splitting 
approaches should not introduce any error when the 
exact solutions of the chemistry and advection operators 


are known. However, when solving numerically the 
advection operator, with an Eulerian advection scheme, 
undesired numerical diffusion will cause the numerical 
advection operator to not commute with the chemistry 
operator (since nonlinear chemical operators do not 
commute with diffusion, as shown in |Hundsdorfer| 
|2003| )) thus signalling the emergence of 
operator splitting errors in the numerical solution of 
equations ([12])-([14]). 

This one-dimensional problem is relevant in realistic 
global 3D simulations since the transport operator is 
solved utilizing Eulerian numerical schemes, and thus, 
giving rise to undesired numerical diffusion that will 
not commute with the time-integration of the chemistry 
operator. Moreover, in regions in the atmosphere where 
the flow is near (2D) divergent-free (due to a well 
stratified atmosphere) and during the night (or day) so 
that chemistry is independent of space, chemistry and 
transport operators may commute locally in space and 
time as in the ID prototype. 

In more complicated circumstances, for example in 
regions of space close to the terminator line (where the 
day and night boundary is), and in Equatorial regions 
where convection makes the atmosphere be far from 
divergent-free conditions, operator splitting errors 
can be expected to be larger since the advection and 
chemistry operators will not commute. 


and Verwer 


3.1. Analytic steady-state solution 

When the chemistry is fast with respect to trans¬ 
port processes, an exact expression can be found for 
the steady-state solution of system ([T2|)-([T4|). For ex¬ 
ample, by choosing k\ = 1000 and £2 = 2000, 

as in ( [Sportisse] |2000| , and introducing the non-stiff 
combined-chemistry operator^ = (NO) 0 3 - 2 NO 2 , 
we can represent a stiff (fast) chemistry operator as the 
quotient xl e f° r a small parameter e. Equations (T2][T4] ) 
can be re-written, as suggested in equation ([9]), as: 


d NO 
dt 


+ u 


d NO 
dx 


d N0 2 
dt 


+ u 


d N0 2 
dx 


d0 3 d0 3 

-h U - 

dt dx 


X 
e ’ 
X 
e ’ 

6 


(15) 

(16) 
(17) 


Here e represents the stiffness of the system and is 
given by the ratio between the slow advection scales 
and the fast chemistry time scales. For example if 
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u ~ 0 ( 1 ) and k t ~ 10 3 , then £ ~ 10 3 . 


The expression of the steady-state solution of sys¬ 
tem is found by introducing the lumped species NO x = 
NO + N0 2 and O x = 0 3 + A0 2 (References, sportisse 
2000) in order to re-write equations (T5j)-(fT7j) as: 


d NO x 
dt 
dO 
dt 
d_ 0 3 

dt 


+ w 


<9 AO 


+ w 


dO 


+ w 


dx 
d O 3 
dx 


X 


X 


= 0, 

(18) 

= 0 , 

(19) 


( 20 ) 


e 


In this new form, and denoting D/Dt = d/dt + u d/dx , 
it can be seen that the lumped species NO x and O x are 
conserved in time, since 


DNO x 

Dt 


= 0 


and 


DNO x 

Dt 


= 0 . 


3.2. Numerical experiments 

We chose to solve equations (|T2])-([T4|) to simulate 
the fate of an instantaneous release containing the three 
chemicals over a 360 km one-dimensional region. The 
constant flow velocity was chosen to resemble realis¬ 
tic atmospheric values of u = 10 m/s. We prescribed a 
computational spatial domain, v € [0, L] for L = 3000 
km, so that the plume would stay within the domain for 
the whole simulation time, t e [0, T] for T — 10 hours, 
and in order to not introduce any errors due to bound¬ 
ary conditions in the numerical advection operator. The 
values of k\ - 1000 and k 2 - 2000 were chosen for 
the stiff chemistry operator. The effective stiffness of 
the chemistry with respect to the transport is 0 ( 10 -2 ) 
since u ~ 0(10). The initial conditions are given by 
NO(x, 0) = NO 2 (x,0) = 0 2 (x, 0) = p(x), where 


p(x) = 


1 if v e [720,1080] 
0 elsewhere. 


As a consequence, for regions where the three species 
are initially present, the exact asymptotic value of the 
concentration of all species, NO\ NO>\, and O 3 , can 
be found explicitly as a function of the initial concen¬ 
tration of the lumped species. This is achieved in two 
steps. First, by expressing the values of the steady state 
concentrations, NO ^ and NO\, as a function of the con¬ 
served lumped species as: 


NO X (0) = NO'+NOl and 0,(0) = 0\+NO\, (21) 


and substituting them in equation (20). The system 
reaches a chemical steady state when x = (NO) O 3 - 
2 N0 2 = 0, or equivalently when 


[NO x { 0) - [0,(0) - 0\W\ ~ 2[0,(0) - Oj] = 0, (22) 


which is a second order equation for the steady state of 
O 3 with solutions given by 



1 (2 + JV0,(0) - 0,(0)) (23) 

\ V(2 + NOM - O,( 0)) 2 + 80,(0) 


And second, the values of NO\ and NO^ can be found 
by substituting the (physically relevant) positive solu¬ 
tion of ( [23] ) in equations $2A) . For time scales r such 
that r » \/k (for k = min(k\,k 2 )), the system will have 
reached chemical steady-state and from then on, equa¬ 
tions ([T8])-([20} (and thus the original system ([12]) -([14])) 
will behave as a transport-only process propagating the 
steady-state concentrations with a constant velocity u. 


In a 10-hour simulation time period, the initial 
release is advected exactly 360 km to the right, 
and the concentrations of all species have reached 
chemical equilibrium. According to expression ( [23] ), 
o\ = NO' = 1.236, and NO\ = 0.764. The exact 
solution at time t = T = 10 hours is explicitly given 
by 0 2 (x,T) = NO(x,T) = 1.236 x p(x - 360) and 
N0 2 (x , T) = 0.764 x p(x - 360). This is our reference 
solution. 


For the numerical simulations, we implemented 
an explicit, second order accurate (in space), one¬ 
dimensional advection-scheme based on the Lax- 
Wendroff method with superbee slope limiters (See 
|LeVeque| ( [2002 ), pp 112 for details), and used for the 
chemistry, the built-in implicit stiff-ODE integrator 
ode23 from Matlab. In order to minimize contributions 
to the numerical error, to the first term in inequality 
from both the advection scheme and chemistry 
integrator, we utilized a very small internal advection 
time step, A t T = 90 seconds, and set the convergence 
relative-tolerance parameter to 10 -3 in the routine 
ode23 (it adaptively chooses a small internal time step 
in order to meet the prescribed 0 . 1 % error convergence 
criterion). 


We solved equations using multiple first 

order Godunov operator splitting approaches (where 
transport and chemistry were evaluated in different 
orders) for multiple operator-splitting time-steps, 
At = 180, 360, 1800 and 3600 seconds, and for 
multiple grid sizes Ax = 22.5, 45, 90, 180 and 
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O concentration (t= 10 hours) 



Ax (m) xIO 5 



Operator spitting time step (s) 

Figure 1: Behavior of numerical error in the one-dimensional 
transport-chemistry system. The top panel shows the analytical “true” 
and numerical solutions at different grid sizes of the system after a 10- 
hour simulation time. The middle panel shows the errors relative to 
the true solution with different grid sizes and operator splitting ap¬ 
proaches. The bottom panel shows the behavior of the relative er- 7 
rors (RRMS) from the two operator splitting approaches, for fixed 
Ax = 180km and different time steps, when compared to the analytic 
solution. 


360 km (the three largest grid sizes were chosen to 
resemble spatial resolutions of 4° x 5°, 2° x 2.5°, and 
1° x 1.25°, in current 3D global CTMs). The results 
of these numerical simulations and the exact solution 
are plotted in the top plot of Figure [T] The numerical 
solutions corresponding to the multiple operator split¬ 
ting approaches, for a given value of Ax, appear as a 
single curve since their differences were smaller than 
the line-width chosen for the plot. 

The quantification of numerical errors was performed 
using the modified relative root mean square (RRMS), 
commonly used in 3D atmospheric chemistry simula¬ 
tions given by 


dJCi) = 


\ 


iy 

M 4^ 


c A - C B 


(24) 


where Cf and C B are the concentrations of species i 
calculated in simulations A and B , respectively, Q is the 
set of grid-boxes where Cf exceeds a threshold a , and 
M is the number of such grid-boxes. We used a =10 -4 , 
thus neglecting concentrations smaller than ~ 0.01% 
with respect to the original concentration. In our 
one-dimensional experiments, simulation A is the exact 
solution, and simulation B was one of the multiple Go¬ 
dunov operator splitting approaches. The second plot 
of Figure [T] shows the quantity d AB = (1 //) 2/ d AB (Ci) for 
i = 3 species, for the multiple values of At and Ax. In 
this plot, the red triangles represent simulations where 
transport was evaluated last, (x - T), and the green 
dots where chemistry was evaluated last (T - x). This 
plot confirms what is observed in the top plot, i.e., the 
fact that the differences across the multiple operator 
splitting approaches, for a given Ax, are very small 
(< 1 %). 


In the bottom plot of Figure [I] we further show the 
values of the numerical error for the two sequences, 
X ~ T and T - x, for Ax = 180 km, for the multiple 
values of the operator splitting time-steps. We found 
this plot to be representative of the behaviour of the 
numerical error for other values of Ax. Note that while 
the differences across the multiple approaches are 
very small, the interesting mathematical behaviour of 
the numerical error, discussed in section |2.1[ can be 
observed. Indeed, the numerical error of the sequences 
T - x> where the chemistry (the stiff process) is 
evaluated last, produce better numerical results than 
their counter parts x ~ T. Moreover, T - x sequences 
appear to be almost insensitive to the magnitude of the 






















operator splitting time-step (the error even seems to 
grow as At —> 0 as reported in Sportisse, 2000) making 
them a preferred choice, since larger operator splitting 
time steps allow faster computations when exploiting 
the intrinsic parallelizable nature of the chemistry 
operator. The quality of results produced by sequences 
where transport is evaluated last, follows the traditional 
behaviour of linear analysis where the numerical error 
decreases as the operator splitting time decreases. Since 
the magnitude of these first order operator splitting 
errors was so small, we chose to not implement higher 
order operator splitting approaches. 


While the bottom plot of Figure [T] shows a clear 
picture of the magnitude of operator splitting errors 
(< 1 %), we performed transport-only simulations in 
order to verify the magnitude of the numerical errors 
coming from the numerical advection scheme itself. 
The results of these simulations are shown in the top 
plot of Figure [2] Note that while the magnitude of 
the concentration of O 3 in these simulations is exactly 
one (since no chemistry is present), the numerically 
simulated profiles, for the different values of Ax, 
look very similar to those in the top plot of Figure [T] 
Indeed, when computing the modified RRMS error 
associated to these simulations, as shown in the bottom 
plot of Figure [2j the behaviour of the relative errors 
resembles the one observed in the middle plot of Figure 
[I] In short, the numerical errors coming from the 
choice of operator splitting are eclipsed by the largest 
component of the numerical error coming from the 
spatial discretization (second term in inequality ([3])) in 
the numerical advection scheme. 


Having chosen an initial condition in the shape of a 
step function in our experiments, caused our second or¬ 
der numerical advection scheme to behave as a first or¬ 
der scheme. Indeed the numerical error decreases close 
to linearly in our numerical experiments when using 
the L 2 -norm instead of the modified RRMS (plot not 
shown). Estimates for the numerical errors, in the form 
of an effective numerical diffusion, Dk , for ID first or¬ 
der numerical advection schemes place their value at 
Dh ~ uAx , where u is the mean flow velocity and Av 
the grid spacing. In our ID experiments, these numer¬ 
ical diffusion is of the order D ~ 10 6 m 2 /s. Numerical 
diffusion in 3D global models ( |Lin and Rood| ( |1996| ); 
Sant illana ( 2013|); [Rastigeyev et al. ( 2007| ); |Wild and 


Prather (2006 ); |Pisso et al.| ( |2QQ9 )) is estimated to be 


around 10 5 - 10 6 m 2 /s. These 3D estimates place our 
one-dimensional experiments within a relevant range. 


0 3 concentration (t=10 hours) Transport only 



Ax (m) xIO 5 


Figure 2: Behavior of numerical error in the one-dimensional 
transport-only system. The top panel shows the analytical “true” and 
numerical solutions at different grid sizes of the system after a 10-hour 
simulation time. The bottom panel shows the errors relative to the true 
solution with different grid sizes and operator splitting approaches. 
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4. Numerical experiments using GEOS-Chem 


Determining the exact magnitude of numerical errors 
in 3D global CTM simulations in the exact same way 
we did for our ID prototype is not possible. This 
is due to the lack of an analytic expression for the 
solution to equations 0 in realistic circumstances 
(time-dependent winds, time-dependent chemistry rates 
changing throughout the geographic domain due to 
photolysis, time-dependent emissions). In order to 
estimate operator splitting errors in 3D CTMs, we 
can only compare the output of simulations where 
everything is kept the same except for the operator 
splitting sequence and the operator splitting time step. 
This is the strategy we present in this section, which in 
combination with the results from our one-dimensional 
simulations, allowed us to determine upper bounds of 
operator splitting errors in GEOS-Chem. In order to 
further understand additional numerical errors, due to 
the time integration of relevant processes inherently 
solved with operator splitting approaches and not 
present in our ID toy example, we performed multiple 
additional 3D global simulations. In these simulations, 
we gradually included inhomogeneous boundary 
conditions (emission processes) to the time integration 
( Sportisse] ( 2000| ); |Hundsdorfer and Verwer| ( |2Q03| )), 
and vertical processes (convection and dry deposition). 


GEOS-Chem is a state-of-the-art 3D global Eule- 
rian model of tropospheric chemistry driven by as¬ 
similated meteorological observations from the God¬ 
dard Earth Observing System (GEOS) of the NASA 
Global Modeling and Assimilation Ofice (GMAO). 
The model simulates global tropospheric ozone-NOx- 
VOC-aerosol chemistry. The full chemical mech¬ 
anism for the troposphere involves over a hundred 
species and over three hundred reactions. The 
ozone- NO x - HO x - VOC-aerosol chemical mecha¬ 
nism of GEOS-Chem has been described by |Bey et| 
al. (|2Q01|); |Park et al| ( |2004| ) and recently updated 
by M ao et al ( |201Q ). Details of the chemical reac¬ 
tions and rate constants are reported in the chemi¬ 
cal mechanism document (http://acmg.seas.harvard.edu 
/geos/wiki_docs/chemistry /chemistry _updates_v6.pdf). 
In Figures 4-6 the chemical species are arranged in 
the order of their chemical lifetimes in the atmosphere, 
from OH (< 1 second) and NO x (~1 hour), to CO and 
C 2#6 (2-3 months). 

The chemical mass balance equations are integrated 
using a Gear-type solver ( Jacobson] 1995| ). Strato¬ 
spheric chemistry is not explicitly simulated and it 
instead uses the Synoz cross-tropopause ozone flux 


boundary condition of McLinden et al. ( |2000 ). The 
model uses the flux form semi-Langrangian advection 
scheme of |Lin and Rood| (l996| ). We used the GEOS- 
Chem model (v8-02-03) driven by the GEOS-5 data at 
the 4 x 5 horizontal resolution and 47 levels in the ver¬ 
tical. Detailed descriptions of the model are given by 
( |Bey et al.]|2001| ) and ( [Zhang et al||2011| ). In this study, 
we initiate the model simulations on January 1, 2005 
with model fields from a 6-month spin-up run, and fo¬ 
cus on the weekly averaged model results for January 
1-7, 2005. 


4.1. Transport and chemistry 

Our strategy consisted of comparing the instanta¬ 
neous concentration of several chemical species, after 
multiple one-week long, 4° x 5° horizontal resolution, 
GEOS-Chem simulations (version v8-02-02), using two 
versions of the (default) second order Strang operator 
splitting method given by the sequences: 

T(At/2) X (At)T(At/2) and X (At/2)T(At) X (At/2) 

for different values of the operator-splitting time step 
At. These sequences are denoted as T X T and X T X 
respectively in the subsequent paragraphs. We used 
At = 60, 30,10,2 mins. In all these simulations, trans¬ 
port and chemistry were the only active mechanisms, 
all other mechanisms were turned off. The inactive 
mechanisms include: emissions, convective transport, 
deposition, and planetary boundary layer mixing. 
Emissions correspond to inhomogenous boundary 
conditions that are treated numerically as production 
rates distributed in the boundary layer and solved 
together in the chemistry operator. 

We used the modified RRMS ( [24] ) with a threshold 
a - 10 6 molecules cm -3 to quantify the numerical dif¬ 
ferences in our global simulations. Figure [4] shows the 
relative differences between the reference simulation 
X T X with At - 2 mins, and the other operator splitting 
approaches for multiple Af’s. Note that the maximum 
differences across simulations (and species) are of the 
order of ~ 10%. 

Using our one-dimensional prototype and fixing 
Av = 180 km, we compared the results of two operator 
splitting strategies ( T — X &vA X - T) for multiple values 
of At, with the sequence T — X and At = 3600 sec 
set as a reference. The results are displayed in Figure 
[3] Note that while the bottom plot of Figure [I] shows 
that operator splitting (relative) errors are less than 1% 
(when comparing to the analytic solution), the relative 
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differences between simulations using alternative oper¬ 
ator splitting methods may be as large as 10%. This is 
roughly the same magnitude of the differences observed 
between the 3D (transport-chemistry) simulations in 
the top panel of Figure]?] 

Note that we chose the sequence x^X with Af = 2 
mins as the reference simulation for our 3D experi¬ 
ments, instead of the sequence x^X with At = 60 that 
would have been suggested by our ID experiments 
(as in Figure [3j. The reason for this is shown in 
the top panel of Figure [5] where we can see that the 
differences between (transport-chemistry) simulations 
with different operator splitting sequences but with 
the same time step, get smaller as At gets smaller. 
This behaviour would be expected from a converging 
operator splitting method where none of the operators is 
stiff and where the order of evaluation of the operators 
is not relevant. An alternative explanation could be 
that the operator splitting errors are very small and 
what we are observing is the convergence of the 
time-integration of each operator, separately, as At 
gets small. This would suggest that the numerical 
errors of the time-integration of the transport and 
the chemistry contribute significantly to the first term 
(involving At) on the right hand-side of inequality 
and should be comparable, in magnitude, to those 
observed between different operator splitting sequences. 

In order to investigate this, we plotted the differences 
between simulations where the only active mechanism 
was either chemistry or transport, for multiple AFs, 
while keeping all other parameters exactly the same 
as in the previous simulations. The results are plotted 
in Figure [6] These two plots show that indeed the 
numerical errors arising from the time-integration of 
each of the operators separately lead to differences of 
the same magnitude as those observed in the operator 
splitting simulations. We also observe that the differ¬ 
ences get smaller as At decreases suggesting numerical 
convergence. This comparable differences make it hard 
to disentangle a sharp estimate of the operator splitting 
in 3D. 

Note also that in our one dimensional prototype a 
cleaner analysis was achieved since we chose a smaller 
internal time step (A t T = 90 seconds) to integrate the 
(explicit-in-time) transport operator than the operator 
splitting time step (180 seconds< At < 60mins). 
This choice reduced the contribution to the numerical 
errors involving At in inequality Q from the transport 
integration. In order to save computational time in 


GEOS-Chem (and in most CTMs), however, the time 
step of the (explicit-in-time) transport scheme is chosen 
to be equal to the operator splitting time step leading to 
larger numerical errors. 


In our one dimensional prototype, the chemistry 
operator was solved using an adaptive time-integration 
routine with very tight convergence constraints, thus 
reducing numerical errors. The time-integration of the 
chemistry operator in GEOS-Chem uses an adaptive 
time stepping strategy ( Jacobson (1995])) in order to 
meet convergence requirements (absolute and relative 
numerical error tolerances) at every user-defined time 
step. These parameters have been internally set to 
keep simulation times reasonable while maintaining 
acceptable numerical accuracy. We kept these settings 
as they are typically used in global simulations for our 
numerical experiments. Figure [6] shows the differences 
between chemistry-only simulations for different user- 
defined chemistry time-steps. Presumably these errors 
could be decreased by fine tuning the error tolerances 
in the time integration routine appropriately, but this 
approach may increase processing times considerably. 


Despite all of these numerical issues, we highlight 
the fact that we can establish an upper limit of about a 
10% for the magnitude of operator splitting errors based 
on the results of our multiple simulations in 3D. More¬ 
over, we show that differences of the single chemical 
species with largest discrepancies across simulations, 
Isoprene, are not significant in Figures [8] [7] and [9] for 
chemistry only simulations, transport only simulations, 
and different sequences of operator splitting methods, 
respectively. From these plots and the results of our 
one-dimensional prototype, we hypothesize that the 
operator splitting errors may be much smaller than 10%. 


We also highlight the fact that we did not pursue fur¬ 
ther efforts to show that the sequences evaluating the 
chemistry at the end of the time step in 3D compare bet¬ 
ter with observations, since our one-dimensional pro¬ 
totype, as well as multiple studies in global CTMs 


([Rastigeyev et al.|[2007||Prather et al.|[200 8; Sanlil lana 


2013), suggest that the numerical errors associated to 
the transport integration, at current spatial resolutions, 
are significantly larger than those observed in operator 
splitting methods. In addition, uncertainties in emis¬ 
sion fields and deposition mechanisms may pose fur¬ 
ther difficulties in addressing this question. In our one¬ 
dimensional proto-type, subsequent reductions in the 
spatial resolution lead to significant improvements in 
the accuracy of the numerical solution globally (for any 
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operator splitting sequence). Whereas a better choice 
of operator splitting (where chemistry is evaluated last) 
leads to a very modest improvement at a given spatial 
resolution Av. 


4.2. Boundary conditions and vertical processes 


Other important processes in 3D simulations are in¬ 
tegrated in time using operator splitting strategies. As 
noted in [Sportisse| ( |2000| ) and Hundsdorfer jmd Verwer 
( 2003) the time integration of inhomogeneous boundary 
conditions, such as emission processes in global simu¬ 
lations, using operator splitting strategies may lead to 
considerable numerical errors. Additionally, the time 
integration of vertical processes such as convection and 
deposition using operator splitting may also lead to im¬ 
portant numerical errors. 


In order to investigate the magnitude of numerical 
errors due to these processes, we performed additional 
3D simulations that gradually included inhomogeneous 
boundary conditions (emissions) and vertical process. 
In other words, aside from the 3D “transport-chemistry” 
simulations discussed in the previous sections, we per¬ 
formed simulations with (i) “transport, chemistry, and 
emissions” and simulations with (ii) “transport, chem¬ 
istry, emissions, convective transport and deposition”. 
When emissions are included, they are integrated within 
the chemistry solver, using the chemistry time step. 
Convective transport and deposition are solved using the 
standard setting of GEOS-Chem, which integrate these 
two processes (sequentially) during the chemistry time 
step. 


The differences between these two sets of simula¬ 
tions, using the same methodology explained in the pre¬ 
vious section, are plotted in the two lower panels of Fig¬ 
ures [4] and [5] As these Figures show, the additional nu¬ 
merical errors coming from the inclusion of inhomoge- 
nous boundary conditions (emissions) are significant. 
Indeed, the differences between the simulations that in¬ 
clude “transport, chemistry, and emissions” are roughly 
double the magnitude of the differences between the 
simulations that include only “transport-chemistry” for 
different operator splitting strategies. The incorpora¬ 
tion of convective transport and deposition to the sim¬ 
ulations does increase the differences between simula¬ 
tion, mainly when changes in time steps are large, as 
shown in the bottom panel of Figure]?] When time-steps 
are fixed and operator splitting approaches are different, 
these vertical processes do not seem to lead to larger 
differences in the different simulations. 



Figure 3: Behavior of the relative errors (RRMS) of simulations per¬ 
formed with two different operator splitting approaches (T - x an d 
X - T ), fixing Ax = 180 km, for multiples time steps. The reference 
solution is obtained with the sequence T ~x for At = 3600 seconds. 


5. Conclusions and Future work 


We have presented a way to characterize operator 
splitting errors in the context of atmospheric chem¬ 
istry modeling. Our approach numerically extends one¬ 
dimensional linear results to non-linear ID and 3D 
cases. These numerical findings are relevant to global 
atmospheric chemistry modeling. Our findings suggest 
that stiff operators should be evaluated lastly in oper¬ 
ator splitting methodologies. This results is consistent 
with the linear results presented in ( |Sportisse] |2000[ ), 
and previous studies in numerical weather prediction 


(Dubai et al. 2005). Differences of approximately 10% 
across species are found when comparing the outputs of 
global simulations using different operator splitting ap¬ 
proaches, using multiple splitting time steps. This, in 
combination with our one-dimensional results, suggests 
that operator splitting errors do not exceed 10% rela¬ 
tive errors in global simulations. We show also, that in 
current spatial resolutions, the numerical diffusion er¬ 
rors introduced in global atmospheric chemistry mod¬ 
els eclipse errors emerging from operator splitting tech¬ 
niques. 


5.1. Future work 

Future studies should identify whether operator split¬ 
ting strategies that evaluate fast dynamics operators 
lastly in global simulations lead to simulations that 
improve the match between simulations and observa¬ 
tions. Further exploration is also required regard¬ 
ing the effect of different operator splitting strategies 
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Figure 4: Behavior of numerical error in the GEOS-Chem 3-D model 
simulations. Here TCT denotes Transport-Chemistry-Transport, CTC 
denotes Chemistry-Transport-Chemistry, and the numbers denotes op¬ 
erator splitting time steps in minutes. Relative RMS relative to the 
CTC2 model simulation are shown for different chemical species with 
lifetimes ranging from seconds (OH) to months (CO,C 2 H^). Active 
processes in these simulations are as follows: Transport and chem¬ 
istry (top panel); Transport, chemistry and emissions (middle panel); 
Trasnport, chemistry, emissions, convective transport and deposition 
(bottom panel). 


Figure 5: Behavior of numerical error in the GEOS-Chem 3-D model 
simulations. Here TCT denotes Transport-Chemistry-Transport, CTC 
denotes Chemistry-Transport-Chemistry, and the numbers denotes op¬ 
erator splitting time steps in minutes. Relative RMS for different 
operator splitting approaches for fixed time steps: At = 2,30,60 
mins.Active processes in these simulations are as follows: Transport 
and chemistry (top panel); Transport, chemistry and emissions (mid¬ 
dle panel); Trasnport, chemistry, emissions, convective transport and 
deposition (bottom panel). 
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a) Transport only 



Figure 6: Behavior of numerical error in the GEOS-Chem 3-D model 
simulations. Relative RMS for transport-only (top panel) and chem¬ 
istry only (bottom) simulations using different time steps: At = 
2,30,60 mins. 





Figure 7: Comparison of isoprene concentrations using different time 
steps for GEOS-Chem transport-only simulations. Isoprene concen¬ 
trations at the surface level from the model simulation with time step 
of 60 minutes (top-left panel) are compared to the model simulation 
with time step of 2 minutes (top-right panel). Absolute (bottom-left) 
and relative differences (bottom-right) are also shown. 


in the time integration of of the governing equations 
of aerosol dynamics and different choices of bound¬ 
ary layer mixing schemes. Additional “toy-tests” that 
should be explored in order to further understand nu¬ 
merical errors introduced by different operator splitting 
strategies include those discussed in ( [Lauritzenj |2Q14| 
Pudykiewicz, 2006 ). Finally, nuances between operator 
splitting approaches in Eulerian and Semi-Lagrangian 
transport schemes should be more deeply investigated 
( [Pudykiewicz et ak||l997| ). 
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